Load Libraries & ArchR object

library(ArchR)
library(rhdf5)
library(reticulate)
library(tidyverse)
#library(caret)
h5disableFileLocking()

# load ArchR project
proj <- loadArchRProject(path = "Save_Vignette_object_4/")

# for R to find macs2 we will use a conda environment where the latest version 
# of macs2 is installed
#conda_list()[[1]][2] %>% use_condaenv(required = TRUE)

Create Low-Overlapping Aggregates of Cells

If we compute correlation on sparse features, this can lead to noise in the data. Therefore, first low-overlapping aggregates of single cells are created (see Cicero). Aggregates with >80% overlap with an other aggregate are removed.

  • sampling with replacement, k = 100 (default) groups/cell
  • aggregate binary accessibility profiles for cells in each group

Co-accessibility

Correlation in accessiblity between two peaks across many single cells. For example if peak A is accessible in a single cell is peak B also accessible? With this analysis we might learn if a certain enhancer peak is often co-accessible with a certain promoter peak in a particular cell type.

Keep in mind that cell-type specific peaks would be identified as co-accessible within this cell tpye. This strong correlation does not mean that there is an underlying co-accessibility!

Query hits and subject hits are the indices of the two peaks that are correlated. The correlation column contains the numeric correlation value of the accessibility between the two peaks.

# compute peak co-accessibility information and store it in the ArchR project
proj <- addCoAccessibility(proj, 
                           reducedDims = "IterativeLSI",
                           k = 100,
                           corCutOff = 0.75,
                           # max. 80 % overlap between current group and all previous groups to permit
                          # current grup to be added to the group list during k-nearest neighbor calculation
                          overlapCutoff = 80,
                          dimsToUse = 1:30,
                           #number of k-nearest neighbor groupings to test for passing the supplied overlapCutoff
                          knnIteration = 500 )

# to retrieve the co-accessibility information we use getCoaccessibility()
# return dataframe if returnLoops = FALSE
coacc <- getCoAccessibility(
    ArchRProj = proj,
    # if the dimension has a correaltion to the sequencing depth 
    # that is greater than the cutoff it will be excluded
    corCutOff = 0.5,
    resolution = 1,
    returnLoops = FALSE # if TRUE it returns a GRangesObject
)

coacc %>% head %>% knitr::kable()
queryHits subjectHits seqnames correlation Variability1 Variability2 TStat Pval FDR VarQuantile1 VarQuantile2
6 9 chr1 0.5948088 0.0041591 0.0302828 16.51231 0 0 0.4991376 0.9302221
9 6 chr1 0.5948088 0.0302828 0.0041591 16.51231 0 0 0.9302221 0.4991376
9 31 chr1 0.6909559 0.0302828 0.0141892 21.32989 0 0 0.9302221 0.8024554
18 19 chr1 0.5128650 0.0259368 0.0052909 13.33193 0 0 0.9060243 0.5682542
19 18 chr1 0.5128650 0.0052909 0.0259368 13.33193 0 0 0.5682542 0.9060243
21 49 chr1 0.5345478 0.0042255 0.0335266 14.11476 0 0 0.5038296 0.9449599

Metadata:

metadata(coacc)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
##             seqnames              ranges strand
##                <Rle>           <IRanges>  <Rle>
##        Mono     chr1       752506-753006      *
##           B     chr1       762688-763188      *
##         GMP     chr1       773653-774153      *
##           B     chr1       801007-801507      *
##           B     chr1       805041-805541      *
##         ...      ...                 ...    ...
##        Mono     chrX 154840753-154841253      *
##          NK     chrX 154842390-154842890      *
##          NK     chrX 154862033-154862533      *
##   Erythroid     chrX 154912392-154912892      *
##         GMP     chrX 154997007-154997507      *
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
dim(coacc)
## [1] 111132     11

We can also return the co-accessibility information as a loop track. The resolution parameter sets the base-pair resolution of the loops. Resolution = 1 creates loops that connect the center of each peak.

Where is the information contained which peaks are correlated?

coacc <- getCoAccessibility(ArchRProj = proj,
                            corCutOff = 0.5,
                            resolution = 1,
                            returnLoops = TRUE)


coacc [[1]]
## GRanges object with 55566 ranges and 9 metadata columns:
##           seqnames              ranges strand | correlation Variability1
##              <Rle>           <IRanges>  <Rle> |   <numeric>    <numeric>
##       [1]     chr1       845703-856628      * |    0.594809   0.00415906
##       [2]     chr1       856628-940499      * |    0.690956   0.01418923
##       [3]     chr1       894695-895356      * |    0.512865   0.00529094
##       [4]     chr1       901579-999702      * |    0.534548   0.00422549
##       [5]     chr1      948868-1004761      * |    0.506462   0.02580794
##       ...      ...                 ...    ... .         ...          ...
##   [55562]     chrX 153583564-153665944      * |    0.537871  0.000662602
##   [55563]     chrX 153584852-153597227      * |    0.566064  0.001244477
##   [55564]     chrX 153597227-153637613      * |    0.512202  0.063475852
##   [55565]     chrX 153686270-153744595      * |    0.513243  0.012864089
##   [55566]     chrX 153947827-154027599      * |    0.524475  0.013135018
##           Variability2     TStat        Pval         FDR VarQuantile1
##              <numeric> <numeric>   <numeric>   <numeric>    <numeric>
##       [1]    0.0302828   16.5123 3.62193e-49 2.50640e-47     0.499138
##       [2]    0.0302828   21.3299 3.41000e-72 6.81189e-70     0.802455
##       [3]    0.0259368   13.3319 6.84909e-35 2.23707e-33     0.568254
##       [4]    0.0335266   14.1148 2.71023e-38 1.06988e-36     0.503830
##       [5]    0.0189477   13.1076 6.23013e-34 1.92659e-32     0.905221
##       ...          ...       ...         ...         ...          ...
##   [55562]   0.00116998   14.2381 7.75355e-39 3.14586e-37    0.0623195
##   [55563]   0.06347585   15.3237 1.05928e-43 5.56776e-42    0.1644870
##   [55564]   0.02888357   13.3086 8.62737e-35 2.80265e-33    0.9912758
##   [55565]   0.02225960   13.3453 6.00370e-35 1.96741e-33    0.7838576
##   [55566]   0.00519124   13.7465 1.10713e-36 3.99906e-35    0.7879772
##           VarQuantile2     value
##              <numeric> <numeric>
##       [1]     0.930222  0.594809
##       [2]     0.930222  0.690956
##       [3]     0.906024  0.512865
##       [4]     0.944960  0.534548
##       [5]     0.853012  0.506462
##       ...          ...       ...
##   [55562]     0.151581  0.537871
##   [55563]     0.991276  0.566064
##   [55564]     0.923357  0.512202
##   [55565]     0.880343  0.513243
##   [55566]     0.562910  0.524475
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

If we decrease the resolution to 1000 this could help with overlplotting of co-accessibility interactions. The GRanges object will contain fewer total entries, 51,82 instaed of 54,002.

cA <- getCoAccessibility(
    ArchRProj = proj,
    corCutOff = 0.5,
    resolution = 1000,
    returnLoops = TRUE
)

cA[[1]]
## GRanges object with 53116 ranges and 9 metadata columns:
##           seqnames              ranges strand | correlation Variability1
##              <Rle>           <IRanges>  <Rle> |   <numeric>    <numeric>
##       [1]     chr1       845500-856500      * |    0.594809   0.00415906
##       [2]     chr1       856500-940500      * |    0.690956   0.01418923
##       [3]     chr1       894500-895500      * |    0.512865   0.00529094
##       [4]     chr1       901500-999500      * |    0.534548   0.00422549
##       [5]     chr1      948500-1004500      * |    0.506462   0.02580794
##       ...      ...                 ...    ... .         ...          ...
##   [53112]     chrX 153583500-153665500      * |    0.537871  0.000662602
##   [53113]     chrX 153584500-153597500      * |    0.566064  0.001244477
##   [53114]     chrX 153597500-153637500      * |    0.512202  0.063475852
##   [53115]     chrX 153686500-153744500      * |    0.513243  0.012864089
##   [53116]     chrX 153947500-154027500      * |    0.524475  0.013135018
##           Variability2     TStat        Pval         FDR VarQuantile1
##              <numeric> <numeric>   <numeric>   <numeric>    <numeric>
##       [1]    0.0302828   16.5123 3.62193e-49 2.50640e-47     0.499138
##       [2]    0.0302828   21.3299 3.41000e-72 6.81189e-70     0.802455
##       [3]    0.0259368   13.3319 6.84909e-35 2.23707e-33     0.568254
##       [4]    0.0335266   14.1148 2.71023e-38 1.06988e-36     0.503830
##       [5]    0.0189477   13.1076 6.23013e-34 1.92659e-32     0.905221
##       ...          ...       ...         ...         ...          ...
##   [53112]   0.00116998   14.2381 7.75355e-39 3.14586e-37    0.0623195
##   [53113]   0.06347585   15.3237 1.05928e-43 5.56776e-42    0.1644870
##   [53114]   0.02888357   13.3086 8.62737e-35 2.80265e-33    0.9912758
##   [53115]   0.02225960   13.3453 6.00370e-35 1.96741e-33    0.7838576
##   [53116]   0.00519124   13.7465 1.10713e-36 3.99906e-35    0.7879772
##           VarQuantile2     value
##              <numeric> <numeric>
##       [1]     0.930222  0.594809
##       [2]     0.930222  0.690956
##       [3]     0.906024  0.512865
##       [4]     0.944960  0.534548
##       [5]     0.853012  0.506462
##       ...          ...       ...
##   [53112]     0.151581  0.537871
##   [53113]     0.991276  0.566064
##   [53114]     0.923357  0.512202
##   [53115]     0.880343  0.513243
##   [53116]     0.562910  0.524475
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

Plot Browser Tracks of Co-accessibility

We can add the co-accessibility information as loop tracks to our browser track. For this we can use the loops parameter.

# marker genes we are interested in 
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

# plot
p <- plotBrowserTrack(
    ArchRProj = proj, 
    groupBy = "Clusters2", 
    geneSymbol = markerGenes, 
    upstream = 50000,
    downstream = 50000,
    loops = getCoAccessibility(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
##       seqnames              ranges strand |     gene_id      symbol
##          <Rle>           <IRanges>  <Rle> | <character> <character>
##   [1]     chr1 208059883-208084683      - |         947        CD34
##   [2]     chrX   48644982-48652717      + |        2623       GATA1
##   [3]     chr9   36838531-37034476      - |        5079        PAX5
##   [4]    chr11   60223282-60238225      + |         931       MS4A1
##   [5]     chr5 140011313-140013286      - |         929        CD14
##   [6]    chr11 118209789-118213459      - |         915        CD3D
##   [7]     chr2   87011728-87035519      - |         925        CD8A
##   [8]    chr17   45810610-45823485      + |       30009       TBX21
##   [9]     chr5   35856977-35879705      + |        3575        IL7R
##   -------
##   seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

grid::grid.newpage()
grid::grid.draw(p$CD14)

Peak2Gene Linkage

What is the difference between co-accessibility and peak2gene linkage?

Co-accessibility uses only scATAC-seq data looking for correlations in accesssiblity.

Peak2Gene linkage in comparison also uses integrated scRNA-seq data to look for correlations between peak accessibility and gene expression.

Both approaces are solutions to a similar problem, but peak2gene linkage is assumed to provide more biologically relevant regulatory links.

All potential peak-to-gene linkages are identified and significant links (R > 0.45 and FDR < 0.1) are kept (Higher correlation and variance suggest that a link is not random noise.)

proj <- addPeak2GeneLinks(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  useMatrix = "GeneIntegrationMatrix",
  dimsToUse = 1:30,
  k = 100, # number of k-nearest neighbors to use for creating single cell goups for correlation analysis
  maxDist = 250000, # maximum allowable distance in basepairs between two peaks to consider for co-accessiblity
  scaleTo = 10^4 # total insertion coutns from a group of cell is summed across all peaks and normalized to toal depth provided by scaleTo
)
p2g <- getPeak2GeneLinks(
    ArchRProj = proj,
    corCutOff = 0.45,
    resolution = 1,
    returnLoops = FALSE,
    varCutOffATAC = .25, # minimum variance quantile of ATAC peak accessibility
    varCutOffRNA = .25, # minimum variance quantile of the RNA gene expression when selecting links
    FDRCutOff = 1e-04 # maximum numeric peak2gene FDR to return
)

as_data_frame(p2g) %>% arrange(desc(Correlation)) %>% 
  head %>% knitr::kable(caption = "Peak-to-gene linkage")
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.
Peak-to-gene linkage
idxATAC idxRNA Correlation FDR VarQATAC VarQRNA
76211 2257 0.9776989 0 0.8779933 0.6813612
42485 14306 0.9621724 0 0.9966034 0.9910757
23457 10826 0.9599931 0 0.9372381 0.5969572
81836 3162 0.9599560 0 0.9999645 0.9820440
30183 12832 0.9559511 0 0.2635277 0.5640019
7959 1140 0.9556089 0 0.9458961 0.3378851

In the metadata the peakSet with ranges, strand and chromosome are saved, as well as the gene set with ranges, strand, chromosome and gene name.

metadata(p2g)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
##            seqnames              ranges strand
##               <Rle>           <IRanges>  <Rle>
##        [1]     chr1       752506-753006      *
##        [2]     chr1       762688-763188      *
##        [3]     chr1       773653-774153      *
##        [4]     chr1       801007-801507      *
##        [5]     chr1       805041-805541      *
##        ...      ...                 ...    ...
##   [141021]     chrX 154840753-154841253      *
##   [141022]     chrX 154842390-154842890      *
##   [141023]     chrX 154862033-154862533      *
##   [141024]     chrX 154912392-154912892      *
##   [141025]     chrX 154997007-154997507      *
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
## 
## $geneSet
## GRanges object with 18601 ranges and 2 metadata columns:
##           seqnames    ranges strand |      name     idx
##              <Rle> <IRanges>  <Rle> |   <array> <array>
##       [1]     chr1     69091      * |     OR4F5       1
##       [2]     chr1    762902      * | LINC00115       2
##       [3]     chr1    812182      * |    FAM41C       3
##       [4]     chr1    860530      * |    SAMD11       4
##       [5]     chr1    894679      * |     NOC2L       5
##       ...      ...       ...    ... .       ...     ...
##   [18597]     chrX 154299695      * |     BRCC3     708
##   [18598]     chrX 154444701      * |      VBP1     709
##   [18599]     chrX 154493852      * |    RAB39B     710
##   [18600]     chrX 154563986      * |     CLIC2     711
##   [18601]     chrX 154842622      * |     TMLHE     712
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
## 
## $seATAC
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seATAC-Group-KNN.rds"
## 
## $seRNA
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seRNA-Group-KNN.rds"
atac_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seATAC-Group-KNN.rds")
rna_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seRNA-Group-KNN.rds")
atac_knn
## class: RangedSummarizedExperiment 
## dim: 141025 491 
## metadata(1): KNNList
## assays(2): ATAC RawATAC
## rownames(141025): f1 f2 ... f141024 f141025
## rowData names(13): score replicateScoreQuantile ... idx N
## colnames: NULL
## colData names(0):
dim(assays(atac_knn)[[1]])
## [1] 141025    491
#assays(atac_knn)[[1]] %>% head

#assays(rna_knn)[[1]] %>% head

We can also return loops:

p2g_loop <- getPeak2GeneLinks(
    ArchRProj = proj,
    corCutOff = 0.45,
    resolution = 1000, # decrease resolution for plotting
    returnLoops = TRUE
)

p2g_loop[[1]] %>% head
## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames         ranges strand |     value         FDR
##          <Rle>      <IRanges>  <Rle> | <numeric>   <numeric>
##   [1]     chr1  762500-812500      * |  0.450391 7.66477e-25
##   [2]     chr1  762500-895500      * |  0.460307 4.97788e-26
##   [3]     chr1  762500-948500      * |  0.515106 2.34837e-33
##   [4]     chr1  812500-894500      * |  0.451671 5.41109e-25
##   [5]     chr1  812500-895500      * |  0.496871 9.11031e-31
##   [6]     chr1 812500-1004500      * |  0.460707 4.44943e-26
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

p <- plotBrowserTrack(
    ArchRProj = proj, 
    groupBy = "Clusters2", 
    geneSymbol = markerGenes, 
    upstream = 50000,
    downstream = 50000,
    loops = getPeak2GeneLinks(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
##       seqnames              ranges strand |     gene_id      symbol
##          <Rle>           <IRanges>  <Rle> | <character> <character>
##   [1]     chr1 208059883-208084683      - |         947        CD34
##   [2]     chrX   48644982-48652717      + |        2623       GATA1
##   [3]     chr9   36838531-37034476      - |        5079        PAX5
##   [4]    chr11   60223282-60238225      + |         931       MS4A1
##   [5]     chr5 140011313-140013286      - |         929        CD14
##   [6]    chr11 118209789-118213459      - |         915        CD3D
##   [7]     chr2   87011728-87035519      - |         925        CD8A
##   [8]    chr17   45810610-45823485      + |       30009       TBX21
##   [9]     chr5   35856977-35879705      + |        3575        IL7R
##   -------
##   seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

grid::grid.newpage()
grid::grid.draw(p$CD14)

Heatmap of peak2gene linkage

The rows are clustered using k-means clustering based on value passed to parameter k (default = 25)

plotPeak2GeneHeatmap(proj, groupBy = "Clusters2")
## Warning: did not converge in 10 iterations

LS0tCnRpdGxlOiAiSW50ZWdyYXRpdmUgQW5hbHlzaXMiCm91dHB1dDogCiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiA1CiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQogICAgdGhlbWU6IGNvc21vCiAgICBoaWdobGlnaHQ6IHRleHRtYXRlCi0tLQoKPHN0eWxlPgpib2R5IHsKdGV4dC1hbGlnbjoganVzdGlmeX0KPC9zdHlsZT4KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBUUlVFLCBhdXRvZGVwID0gVFJVRSwgCiAgICAgICAgICAgICAgICAgICAgICBjb2xsYXBzZSA9IFRSVUUsIG1lc3NhZ2UgPSBGQUxTRSkKa25pdHI6Om9wdHNfa25pdCRzZXQocm9vdC5kaXIgPSAiL29taWNzL2dyb3Vwcy9PRTA1MzMvaW50ZXJuYWwva2F0aGFyaW5hL3NjRG9SSS9BcmNoUiIpCnNldHdkKCIvb21pY3MvZ3JvdXBzL09FMDUzMy9pbnRlcm5hbC9rYXRoYXJpbmEvc2NEb1JJL0FyY2hSLyIpCgpzZXQuc2VlZCgxKQpgYGAKCiMgTG9hZCBMaWJyYXJpZXMgJiBBcmNoUiBvYmplY3QKCmBgYHtyfQpsaWJyYXJ5KEFyY2hSKQpsaWJyYXJ5KHJoZGY1KQpsaWJyYXJ5KHJldGljdWxhdGUpCmxpYnJhcnkodGlkeXZlcnNlKQojbGlicmFyeShjYXJldCkKaDVkaXNhYmxlRmlsZUxvY2tpbmcoKQoKIyBsb2FkIEFyY2hSIHByb2plY3QKcHJvaiA8LSBsb2FkQXJjaFJQcm9qZWN0KHBhdGggPSAiU2F2ZV9WaWduZXR0ZV9vYmplY3RfNC8iKQoKIyBmb3IgUiB0byBmaW5kIG1hY3MyIHdlIHdpbGwgdXNlIGEgY29uZGEgZW52aXJvbm1lbnQgd2hlcmUgdGhlIGxhdGVzdCB2ZXJzaW9uIAojIG9mIG1hY3MyIGlzIGluc3RhbGxlZAojY29uZGFfbGlzdCgpW1sxXV1bMl0gJT4lIHVzZV9jb25kYWVudihyZXF1aXJlZCA9IFRSVUUpCgpgYGAKCiMgQ3JlYXRlIExvdy1PdmVybGFwcGluZyBBZ2dyZWdhdGVzIG9mIENlbGxzCgpJZiB3ZSBjb21wdXRlIGNvcnJlbGF0aW9uIG9uIHNwYXJzZSBmZWF0dXJlcywgdGhpcyBjYW4gbGVhZCB0byBub2lzZSBpbiB0aGUgCmRhdGEuIFRoZXJlZm9yZSwgZmlyc3QgbG93LW92ZXJsYXBwaW5nIGFnZ3JlZ2F0ZXMgb2Ygc2luZ2xlIGNlbGxzIGFyZSBjcmVhdGVkCihzZWUgQ2ljZXJvKS4gQWdncmVnYXRlcyB3aXRoID44MCUgb3ZlcmxhcCB3aXRoIGFuIG90aGVyIGFnZ3JlZ2F0ZSBhcmUgcmVtb3ZlZC4gCgoqIHNhbXBsaW5nIHdpdGggcmVwbGFjZW1lbnQsIGsgPSAxMDAgKGRlZmF1bHQpIGdyb3Vwcy9jZWxsCiogYWdncmVnYXRlIGJpbmFyeSBhY2Nlc3NpYmlsaXR5IHByb2ZpbGVzIGZvciBjZWxscyBpbiBlYWNoIGdyb3VwCgoKIyBDby1hY2Nlc3NpYmlsaXR5CgpDb3JyZWxhdGlvbiBpbiBhY2Nlc3NpYmxpdHkgYmV0d2VlbiB0d28gcGVha3MgYWNyb3NzIG1hbnkgc2luZ2xlIGNlbGxzLiBGb3IgZXhhbXBsZQppZiBwZWFrIEEgaXMgYWNjZXNzaWJsZSBpbiBhIHNpbmdsZSBjZWxsIGlzIHBlYWsgQiBhbHNvIGFjY2Vzc2libGU/IFdpdGggdGhpcwphbmFseXNpcyB3ZSBtaWdodCBsZWFybiBpZiBhIGNlcnRhaW4gZW5oYW5jZXIgcGVhayBpcyBvZnRlbiBjby1hY2Nlc3NpYmxlIHdpdGgKYSBjZXJ0YWluIHByb21vdGVyIHBlYWsgaW4gYSBwYXJ0aWN1bGFyIGNlbGwgdHlwZS4gCgpLZWVwIGluIG1pbmQgdGhhdCBjZWxsLXR5cGUgc3BlY2lmaWMgcGVha3Mgd291bGQgYmUgaWRlbnRpZmllZCBhcyBjby1hY2Nlc3NpYmxlCndpdGhpbiB0aGlzIGNlbGwgdHB5ZS4gVGhpcyBzdHJvbmcgY29ycmVsYXRpb24gZG9lcyBub3QgbWVhbiB0aGF0IHRoZXJlIGlzIGFuIHVuZGVybHlpbmcgY28tYWNjZXNzaWJpbGl0eSEKClF1ZXJ5IGhpdHMgYW5kIHN1YmplY3QgaGl0cyBhcmUgdGhlIGluZGljZXMgb2YgdGhlIHR3byBwZWFrcyB0aGF0IGFyZSAKY29ycmVsYXRlZC4gVGhlIGNvcnJlbGF0aW9uIGNvbHVtbiBjb250YWlucyB0aGUgbnVtZXJpYyBjb3JyZWxhdGlvbiB2YWx1ZSBvZiB0aGUKYWNjZXNzaWJpbGl0eSBiZXR3ZWVuIHRoZSB0d28gcGVha3MuIAoKYGBge3IsIHJlc3VsdHM9ImFzaXMifQojIGNvbXB1dGUgcGVhayBjby1hY2Nlc3NpYmlsaXR5IGluZm9ybWF0aW9uIGFuZCBzdG9yZSBpdCBpbiB0aGUgQXJjaFIgcHJvamVjdApwcm9qIDwtIGFkZENvQWNjZXNzaWJpbGl0eShwcm9qLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgcmVkdWNlZERpbXMgPSAiSXRlcmF0aXZlTFNJIiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgayA9IDEwMCwKICAgICAgICAgICAgICAgICAgICAgICAgICAgY29yQ3V0T2ZmID0gMC43NSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIyBtYXguIDgwICUgb3ZlcmxhcCBiZXR3ZWVuIGN1cnJlbnQgZ3JvdXAgYW5kIGFsbCBwcmV2aW91cyBncm91cHMgdG8gcGVybWl0CiAgICAgICAgICAgICAgICAgICAgICAgICAgIyBjdXJyZW50IGdydXAgdG8gYmUgYWRkZWQgdG8gdGhlIGdyb3VwIGxpc3QgZHVyaW5nIGstbmVhcmVzdCBuZWlnaGJvciBjYWxjdWxhdGlvbgogICAgICAgICAgICAgICAgICAgICAgICAgIG92ZXJsYXBDdXRvZmYgPSA4MCwKICAgICAgICAgICAgICAgICAgICAgICAgICBkaW1zVG9Vc2UgPSAxOjMwLAogICAgICAgICAgICAgICAgICAgICAgICAgICAjbnVtYmVyIG9mIGstbmVhcmVzdCBuZWlnaGJvciBncm91cGluZ3MgdG8gdGVzdCBmb3IgcGFzc2luZyB0aGUgc3VwcGxpZWQgb3ZlcmxhcEN1dG9mZgogICAgICAgICAgICAgICAgICAgICAgICAgIGtubkl0ZXJhdGlvbiA9IDUwMCApCgojIHRvIHJldHJpZXZlIHRoZSBjby1hY2Nlc3NpYmlsaXR5IGluZm9ybWF0aW9uIHdlIHVzZSBnZXRDb2FjY2Vzc2liaWxpdHkoKQojIHJldHVybiBkYXRhZnJhbWUgaWYgcmV0dXJuTG9vcHMgPSBGQUxTRQpjb2FjYyA8LSBnZXRDb0FjY2Vzc2liaWxpdHkoCiAgICBBcmNoUlByb2ogPSBwcm9qLAogICAgIyBpZiB0aGUgZGltZW5zaW9uIGhhcyBhIGNvcnJlYWx0aW9uIHRvIHRoZSBzZXF1ZW5jaW5nIGRlcHRoIAogICAgIyB0aGF0IGlzIGdyZWF0ZXIgdGhhbiB0aGUgY3V0b2ZmIGl0IHdpbGwgYmUgZXhjbHVkZWQKICAgIGNvckN1dE9mZiA9IDAuNSwKICAgIHJlc29sdXRpb24gPSAxLAogICAgcmV0dXJuTG9vcHMgPSBGQUxTRSAjIGlmIFRSVUUgaXQgcmV0dXJucyBhIEdSYW5nZXNPYmplY3QKKQoKY29hY2MgJT4lIGhlYWQgJT4lIGtuaXRyOjprYWJsZSgpCmBgYAoKTWV0YWRhdGE6CgpgYGB7cn0KbWV0YWRhdGEoY29hY2MpCmRpbShjb2FjYykKYGBgCgoKV2UgY2FuIGFsc28gcmV0dXJuIHRoZSBjby1hY2Nlc3NpYmlsaXR5IGluZm9ybWF0aW9uIGFzIGEgbG9vcCB0cmFjay4gVGhlIApyZXNvbHV0aW9uIHBhcmFtZXRlciBzZXRzIHRoZSBiYXNlLXBhaXIgcmVzb2x1dGlvbiBvZiB0aGUgbG9vcHMuIFJlc29sdXRpb24gPSAxCmNyZWF0ZXMgbG9vcHMgdGhhdCBjb25uZWN0IHRoZSBjZW50ZXIgb2YgZWFjaCBwZWFrLiAKCldoZXJlIGlzIHRoZSBpbmZvcm1hdGlvbiBjb250YWluZWQgd2hpY2ggcGVha3MgYXJlIGNvcnJlbGF0ZWQ/CgoKYGBge3J9CmNvYWNjIDwtIGdldENvQWNjZXNzaWJpbGl0eShBcmNoUlByb2ogPSBwcm9qLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgY29yQ3V0T2ZmID0gMC41LAogICAgICAgICAgICAgICAgICAgICAgICAgICAgcmVzb2x1dGlvbiA9IDEsCiAgICAgICAgICAgICAgICAgICAgICAgICAgICByZXR1cm5Mb29wcyA9IFRSVUUpCgoKY29hY2MgW1sxXV0KYGBgCgpJZiB3ZSBkZWNyZWFzZSB0aGUgcmVzb2x1dGlvbiB0byAxMDAwIHRoaXMgY291bGQgaGVscCB3aXRoIG92ZXJscGxvdHRpbmcgb2YgCmNvLWFjY2Vzc2liaWxpdHkgaW50ZXJhY3Rpb25zLiBUaGUgR1JhbmdlcyBvYmplY3Qgd2lsbCBjb250YWluIGZld2VyIHRvdGFsCmVudHJpZXMsIDUxLDgyIGluc3RhZWQgb2YgNTQsMDAyLgoKYGBge3J9CmNBIDwtIGdldENvQWNjZXNzaWJpbGl0eSgKICAgIEFyY2hSUHJvaiA9IHByb2osCiAgICBjb3JDdXRPZmYgPSAwLjUsCiAgICByZXNvbHV0aW9uID0gMTAwMCwKICAgIHJldHVybkxvb3BzID0gVFJVRQopCgpjQVtbMV1dCmBgYAoKCiMjIFBsb3QgQnJvd3NlciBUcmFja3Mgb2YgQ28tYWNjZXNzaWJpbGl0eQoKV2UgY2FuIGFkZCB0aGUgY28tYWNjZXNzaWJpbGl0eSBpbmZvcm1hdGlvbiBhcyBsb29wIHRyYWNrcyB0byBvdXIgYnJvd3NlciAKdHJhY2suIEZvciB0aGlzIHdlIGNhbiB1c2UgdGhlIGBsb29wc2AgcGFyYW1ldGVyLiAKCmBgYHtyfQojIG1hcmtlciBnZW5lcyB3ZSBhcmUgaW50ZXJlc3RlZCBpbiAKbWFya2VyR2VuZXMgIDwtIGMoCiAgICAiQ0QzNCIsICNFYXJseSBQcm9nZW5pdG9yCiAgICAiR0FUQTEiLCAjRXJ5dGhyb2lkCiAgICAiUEFYNSIsICJNUzRBMSIsICNCLUNlbGwgVHJhamVjdG9yeQogICAgIkNEMTQiLCAjTW9ub2N5dGVzCiAgICAiQ0QzRCIsICJDRDhBIiwgIlRCWDIxIiwgIklMN1IiICNUQ2VsbHMKICApCgojIHBsb3QKcCA8LSBwbG90QnJvd3NlclRyYWNrKAogICAgQXJjaFJQcm9qID0gcHJvaiwgCiAgICBncm91cEJ5ID0gIkNsdXN0ZXJzMiIsIAogICAgZ2VuZVN5bWJvbCA9IG1hcmtlckdlbmVzLCAKICAgIHVwc3RyZWFtID0gNTAwMDAsCiAgICBkb3duc3RyZWFtID0gNTAwMDAsCiAgICBsb29wcyA9IGdldENvQWNjZXNzaWJpbGl0eShwcm9qKQopCgpncmlkOjpncmlkLm5ld3BhZ2UoKQpncmlkOjpncmlkLmRyYXcocCRDRDE0KQpgYGAKCgojIFBlYWsyR2VuZSBMaW5rYWdlCgpXaGF0IGlzIHRoZSBkaWZmZXJlbmNlIGJldHdlZW4gY28tYWNjZXNzaWJpbGl0eSBhbmQgcGVhazJnZW5lIGxpbmthZ2U/IAoKQ28tYWNjZXNzaWJpbGl0eSB1c2VzIG9ubHkgc2NBVEFDLXNlcSBkYXRhIGxvb2tpbmcgZm9yIGNvcnJlbGF0aW9ucyBpbiAKYWNjZXNzc2libGl0eS4gCgpQZWFrMkdlbmUgbGlua2FnZSBpbiBjb21wYXJpc29uIGFsc28gdXNlcyBpbnRlZ3JhdGVkIHNjUk5BLXNlcSBkYXRhIHRvIGxvb2sKZm9yIGNvcnJlbGF0aW9ucyBiZXR3ZWVuIHBlYWsgYWNjZXNzaWJpbGl0eSBhbmQgZ2VuZSBleHByZXNzaW9uLiAKCkJvdGggYXBwcm9hY2VzIGFyZSBzb2x1dGlvbnMgdG8gYSBzaW1pbGFyIHByb2JsZW0sIGJ1dCBwZWFrMmdlbmUgbGlua2FnZSBpcyAKYXNzdW1lZCB0byBwcm92aWRlIG1vcmUgYmlvbG9naWNhbGx5IHJlbGV2YW50IHJlZ3VsYXRvcnkgbGlua3MuCgoKQWxsIHBvdGVudGlhbCBwZWFrLXRvLWdlbmUgbGlua2FnZXMgYXJlIGlkZW50aWZpZWQgYW5kIHNpZ25pZmljYW50IGxpbmtzIAooUiA+IDAuNDUgYW5kIEZEUiA8IDAuMSkgYXJlIGtlcHQgKEhpZ2hlciBjb3JyZWxhdGlvbiBhbmQgdmFyaWFuY2Ugc3VnZ2VzdCB0aGF0IAphIGxpbmsgaXMgbm90IHJhbmRvbSBub2lzZS4pCgpgYGB7cn0KcHJvaiA8LSBhZGRQZWFrMkdlbmVMaW5rcygKICBBcmNoUlByb2ogPSBwcm9qLAogIHJlZHVjZWREaW1zID0gIkl0ZXJhdGl2ZUxTSSIsCiAgdXNlTWF0cml4ID0gIkdlbmVJbnRlZ3JhdGlvbk1hdHJpeCIsCiAgZGltc1RvVXNlID0gMTozMCwKICBrID0gMTAwLCAjIG51bWJlciBvZiBrLW5lYXJlc3QgbmVpZ2hib3JzIHRvIHVzZSBmb3IgY3JlYXRpbmcgc2luZ2xlIGNlbGwgZ291cHMgZm9yIGNvcnJlbGF0aW9uIGFuYWx5c2lzCiAgbWF4RGlzdCA9IDI1MDAwMCwgIyBtYXhpbXVtIGFsbG93YWJsZSBkaXN0YW5jZSBpbiBiYXNlcGFpcnMgYmV0d2VlbiB0d28gcGVha3MgdG8gY29uc2lkZXIgZm9yIGNvLWFjY2Vzc2libGl0eQogIHNjYWxlVG8gPSAxMF40ICMgdG90YWwgaW5zZXJ0aW9uIGNvdXRucyBmcm9tIGEgZ3JvdXAgb2YgY2VsbCBpcyBzdW1tZWQgYWNyb3NzIGFsbCBwZWFrcyBhbmQgbm9ybWFsaXplZCB0byB0b2FsIGRlcHRoIHByb3ZpZGVkIGJ5IHNjYWxlVG8KKQpgYGAKCmBgYHtyLCByZXN1bHRzID0gImFzaXMifQpwMmcgPC0gZ2V0UGVhazJHZW5lTGlua3MoCiAgICBBcmNoUlByb2ogPSBwcm9qLAogICAgY29yQ3V0T2ZmID0gMC40NSwKICAgIHJlc29sdXRpb24gPSAxLAogICAgcmV0dXJuTG9vcHMgPSBGQUxTRSwKICAgIHZhckN1dE9mZkFUQUMgPSAuMjUsICMgbWluaW11bSB2YXJpYW5jZSBxdWFudGlsZSBvZiBBVEFDIHBlYWsgYWNjZXNzaWJpbGl0eQogICAgdmFyQ3V0T2ZmUk5BID0gLjI1LCAjIG1pbmltdW0gdmFyaWFuY2UgcXVhbnRpbGUgb2YgdGhlIFJOQSBnZW5lIGV4cHJlc3Npb24gd2hlbiBzZWxlY3RpbmcgbGlua3MKICAgIEZEUkN1dE9mZiA9IDFlLTA0ICMgbWF4aW11bSBudW1lcmljIHBlYWsyZ2VuZSBGRFIgdG8gcmV0dXJuCikKCmFzX2RhdGFfZnJhbWUocDJnKSAlPiUgYXJyYW5nZShkZXNjKENvcnJlbGF0aW9uKSkgJT4lIAogIGhlYWQgJT4lIGtuaXRyOjprYWJsZShjYXB0aW9uID0gIlBlYWstdG8tZ2VuZSBsaW5rYWdlIikKYGBgCkluIHRoZSBtZXRhZGF0YSB0aGUgcGVha1NldCB3aXRoIHJhbmdlcywgc3RyYW5kIGFuZCBjaHJvbW9zb21lIGFyZSBzYXZlZCwKYXMgd2VsbCBhcyB0aGUgZ2VuZSBzZXQgd2l0aCByYW5nZXMsIHN0cmFuZCwgY2hyb21vc29tZSBhbmQgZ2VuZSBuYW1lLgoKYGBge3J9Cm1ldGFkYXRhKHAyZykKYGBgCgpgYGB7cn0KYXRhY19rbm4gPC0gcmVhZFJEUygiU2F2ZV9WaWduZXR0ZV9vYmplY3RfNS9QZWFrMkdlbmVMaW5rcy9zZUFUQUMtR3JvdXAtS05OLnJkcyIpCnJuYV9rbm4gPC0gcmVhZFJEUygiU2F2ZV9WaWduZXR0ZV9vYmplY3RfNS9QZWFrMkdlbmVMaW5rcy9zZVJOQS1Hcm91cC1LTk4ucmRzIikKYXRhY19rbm4KYGBgCgpgYGB7cn0KZGltKGFzc2F5cyhhdGFjX2tubilbWzFdXSkKI2Fzc2F5cyhhdGFjX2tubilbWzFdXSAlPiUgaGVhZAoKI2Fzc2F5cyhybmFfa25uKVtbMV1dICU+JSBoZWFkCmBgYApXZSBjYW4gYWxzbyByZXR1cm4gbG9vcHM6CgpgYGB7cn0KcDJnX2xvb3AgPC0gZ2V0UGVhazJHZW5lTGlua3MoCiAgICBBcmNoUlByb2ogPSBwcm9qLAogICAgY29yQ3V0T2ZmID0gMC40NSwKICAgIHJlc29sdXRpb24gPSAxMDAwLCAjIGRlY3JlYXNlIHJlc29sdXRpb24gZm9yIHBsb3R0aW5nCiAgICByZXR1cm5Mb29wcyA9IFRSVUUKKQoKcDJnX2xvb3BbWzFdXSAlPiUgaGVhZApgYGAKCmBgYHtyfQptYXJrZXJHZW5lcyAgPC0gYygKICAgICJDRDM0IiwgI0Vhcmx5IFByb2dlbml0b3IKICAgICJHQVRBMSIsICNFcnl0aHJvaWQKICAgICJQQVg1IiwgIk1TNEExIiwgI0ItQ2VsbCBUcmFqZWN0b3J5CiAgICAiQ0QxNCIsICNNb25vY3l0ZXMKICAgICJDRDNEIiwgIkNEOEEiLCAiVEJYMjEiLCAiSUw3UiIgI1RDZWxscwogICkKCnAgPC0gcGxvdEJyb3dzZXJUcmFjaygKICAgIEFyY2hSUHJvaiA9IHByb2osIAogICAgZ3JvdXBCeSA9ICJDbHVzdGVyczIiLCAKICAgIGdlbmVTeW1ib2wgPSBtYXJrZXJHZW5lcywgCiAgICB1cHN0cmVhbSA9IDUwMDAwLAogICAgZG93bnN0cmVhbSA9IDUwMDAwLAogICAgbG9vcHMgPSBnZXRQZWFrMkdlbmVMaW5rcyhwcm9qKQopCgpncmlkOjpncmlkLm5ld3BhZ2UoKQpncmlkOjpncmlkLmRyYXcocCRDRDE0KQpgYGAKCiMjIEhlYXRtYXAgb2YgcGVhazJnZW5lIGxpbmthZ2UKClRoZSByb3dzIGFyZSBjbHVzdGVyZWQgdXNpbmcgay1tZWFucyBjbHVzdGVyaW5nIGJhc2VkIG9uIHZhbHVlIHBhc3NlZCB0byAKcGFyYW1ldGVyIGsgKGRlZmF1bHQgPSAyNSkKCmBgYHtyLCBmaWcuaGVpZ2h0PTgsIGZpZy53aWR0aD0xMH0KcGxvdFBlYWsyR2VuZUhlYXRtYXAocHJvaiwgZ3JvdXBCeSA9ICJDbHVzdGVyczIiKQpgYGAKCg==